!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief function that build the poisson section of the input
!> \par History
!>      03.2006 fusing of poisson_dft and poisson_mm
!> \author fawzi
! **************************************************************************************************
MODULE input_cp2k_poisson
   USE bibliography,                    ONLY: &
        Aguado2003, BaniHashemian2016, Blochl1995, Darden1993, Essmann1995, Ewald1921, &
        Genovese2006, Genovese2007, Laino2008, Martyna1999, Toukmaji1996
   USE cell_types,                      ONLY: use_perd_none,&
                                              use_perd_x,&
                                              use_perd_xy,&
                                              use_perd_xyz,&
                                              use_perd_xz,&
                                              use_perd_y,&
                                              use_perd_yz,&
                                              use_perd_z
   USE cp_output_handling,              ONLY: add_last_numeric,&
                                              cp_print_key_section_create,&
                                              low_print_level,&
                                              medium_print_level
   USE cp_units,                        ONLY: cp_unit_to_cp2k
   USE dct,                             ONLY: neumannX,&
                                              neumannXY,&
                                              neumannXYZ,&
                                              neumannXZ,&
                                              neumannY,&
                                              neumannYZ,&
                                              neumannZ
   USE dielectric_types,                ONLY: &
        derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft, derivative_fft_use_deps, &
        derivative_fft_use_drho, rho_dependent, spatially_dependent, spatially_rho_dependent
   USE dirichlet_bc_types,              ONLY: CIRCUMSCRIBED,&
                                              INSCRIBED,&
                                              x_axis,&
                                              xy_plane,&
                                              xz_plane,&
                                              y_axis,&
                                              yz_plane,&
                                              z_axis
   USE input_constants,                 ONLY: do_fist_pol_cg,&
                                              do_fist_pol_none,&
                                              do_fist_pol_sc
   USE input_cp2k_rsgrid,               ONLY: create_rsgrid_section
   USE input_keyword_types,             ONLY: keyword_create,&
                                              keyword_release,&
                                              keyword_type
   USE input_section_types,             ONLY: section_add_keyword,&
                                              section_add_subsection,&
                                              section_create,&
                                              section_release,&
                                              section_type
   USE input_val_types,                 ONLY: enum_t,&
                                              integer_t,&
                                              real_t
   USE kinds,                           ONLY: dp
   USE multipole_types,                 ONLY: do_multipole_charge,&
                                              do_multipole_dipole,&
                                              do_multipole_none,&
                                              do_multipole_quadrupole
   USE ps_implicit_types,               ONLY: MIXED_BC,&
                                              MIXED_PERIODIC_BC,&
                                              NEUMANN_BC,&
                                              PERIODIC_BC
   USE pw_poisson_types,                ONLY: &
        do_ewald_ewald, do_ewald_none, do_ewald_pme, do_ewald_spme, pw_poisson_analytic, &
        pw_poisson_implicit, pw_poisson_mt, pw_poisson_multipole, pw_poisson_periodic, &
        pw_poisson_wavelet
   USE pw_spline_utils,                 ONLY: no_precond,&
                                              precond_spl3_1,&
                                              precond_spl3_2,&
                                              precond_spl3_3,&
                                              precond_spl3_aint,&
                                              precond_spl3_aint2
   USE string_utilities,                ONLY: s2a
#include "./base/base_uses.f90"

   IMPLICIT NONE
   PRIVATE

   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_cp2k_poisson'

   PUBLIC :: create_poisson_section, &
             create_gspace_interp_section, &
             create_ewald_section
!***
CONTAINS

! **************************************************************************************************
!> \brief Creates the Poisson section
!> \param section the section to create
!> \author teo
! **************************************************************************************************
   SUBROUTINE create_poisson_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: subsection

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="poisson", &
                          description="Sets up the poisson resolutor.", &
                          n_keywords=1, n_subsections=0, repeats=.FALSE.)

      NULLIFY (keyword, subsection)
      CALL keyword_create(keyword, __LOCATION__, name="POISSON_SOLVER", &
                          variants=(/"POISSON", "PSOLVER"/), &
                          description="Specify which kind of solver to use to solve the Poisson equation.", &
                          usage="POISSON_SOLVER char", &
                          enum_c_vals=s2a("PERIODIC", "ANALYTIC", "MT", "MULTIPOLE", "WAVELET", "IMPLICIT"), &
                          enum_i_vals=(/pw_poisson_periodic, pw_poisson_analytic, pw_poisson_mt, pw_poisson_multipole, &
                                        pw_poisson_wavelet, pw_poisson_implicit/), &
                          enum_desc=s2a("PERIODIC is only available for fully (3D) periodic systems.", &
                                        "ANALYTIC is available for 0D, 1D and 2D periodic solutions using analytical green "// &
                                        "functions in the g space (slow convergence).", &
                                        "MT (Martyna Tuckermann) decoupling that interacts only with the nearest "// &
                                        "neighbor. Beware results are completely wrong if the cell is smaller than twice the "// &
                                        "cluster size (with electronic density). Available for 0D and 2D systems.", &
                                        "MULTIPOLE uses a scheme that fits the total charge with one gaussian per atom. "// &
                                        "Available only for cluster (0D) systems.", &
                                        "WAVELET allows for 0D, 2D (but only PERIODIC XZ) and 3D systems. It does not "// &
                                        "require very large unit cells, only that the density goes to zero on the faces of "// &
                                        "the cell. The use of PREFERRED_FFT_LIBRARY FFTSG is required.", &
                                        "IMPLICIT allows for 0D, 1D, 2D and 3D systems."), &
                          citations=(/Blochl1995, Martyna1999, Genovese2006, Genovese2007/), &
                          default_i_val=pw_poisson_periodic)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PERIODIC", &
                          description="Specify the directions in which PBC apply. Important notice,"// &
                          " this only applies to the electrostatics."// &
                          " See the CELL section to specify the periodicity used for e.g. the pair lists."// &
                          " Typically the settings should be the same.", &
                          usage="PERIODIC (x|y|z|xy|xz|yz|xyz|none)", &
                          enum_c_vals=s2a("x", "y", "z", "xy", "xz", "yz", "xyz", "none"), &
                          enum_i_vals=(/use_perd_x, use_perd_y, use_perd_z, &
                                        use_perd_xy, use_perd_xz, use_perd_yz, &
                                        use_perd_xyz, use_perd_none/), &
                          default_i_val=use_perd_xyz)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL create_mt_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_wavelet_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_multipole_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_ewald_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_implicit_ps_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)
   END SUBROUTINE create_poisson_section

! **************************************************************************************************
!> \brief Section to set-up parameters for decoupling using the Bloechl scheme
!> \param section the section to create
!> \author teo
! **************************************************************************************************
   SUBROUTINE create_multipole_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: subsection

      CPASSERT(.NOT. ASSOCIATED(section))

      CALL section_create(section, __LOCATION__, name="MULTIPOLE", &
                          description="This section is used to set up the decoupling of QM periodic images with "// &
                          "the use of density derived atomic point charges.", &
                          n_keywords=1, n_subsections=0, repeats=.FALSE.)

      NULLIFY (keyword, subsection)
      CALL keyword_create(keyword, __LOCATION__, name="RCUT", &
                          description="Real space cutoff for the Ewald sum.", &
                          usage="RCUT {real}", n_var=1, type_of_var=real_t, &
                          unit_str="angstrom")
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="EWALD_PRECISION", &
                          description="Precision achieved in the Ewald sum.", &
                          usage="EWALD_PRECISION {real}", n_var=1, type_of_var=real_t, &
                          unit_str="hartree", default_r_val=1.0E-6_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="ANALYTICAL_GTERM", &
                          description="Evaluates the Gterm in the Ewald Scheme analytically instead of using Splines.", &
                          usage="ANALYTICAL_GTERM <LOGICAL>", &
                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="NGRIDS", &
                          description="Specifies the number of grid points used for the Interpolation of the G-space term", &
                          usage="NGRIDS <integer> <integer> <integer> ", n_var=3, default_i_vals=(/50, 50, 50/))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL create_gspace_interp_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL cp_print_key_section_create(subsection, __LOCATION__, "check_spline", &
                                       description="Controls the checking of the G-space term Spline Interpolation.", &
                                       print_level=medium_print_level, filename="GSpace-SplInterp")
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL cp_print_key_section_create(subsection, __LOCATION__, "program_run_info", &
                                       description="Controls the printing of basic information during the run", &
                                       print_level=low_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

   END SUBROUTINE create_multipole_section

! **************************************************************************************************
!> \brief Creates the Martyna-Tuckerman section
!> \param section the section to create
!> \author teo
! **************************************************************************************************
   SUBROUTINE create_mt_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="mt", &
                          description="Sets up parameters of  Martyna-Tuckerman poisson solver. "// &
                          "Note that exact results are only guaranteed if the unit cell is "// &
                          "twice as large as charge density (and serious artefacts can result "// &
                          "if the cell is much smaller).", &
                          n_keywords=1, n_subsections=0, repeats=.FALSE., &
                          citations=(/Martyna1999/))

      NULLIFY (keyword)

      CALL keyword_create(keyword, __LOCATION__, name="ALPHA", &
                          description="Convergence parameter ALPHA*RMIN. Default value 7.0", &
                          usage="ALPHA real", &
                          n_var=1, default_r_val=7.0_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="REL_CUTOFF", &
                          description="Specify the multiplicative factor for the CUTOFF keyword in MULTI_GRID"// &
                          " section. The result gives the cutoff at which the 1/r non-periodic FFT3D is evaluated."// &
                          " Default is 2.0", &
                          usage="REL_CUTOFF real", &
                          n_var=1, default_r_val=2.0_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_mt_section

! **************************************************************************************************
!> \brief ...
!> \param section will contain the ewald section
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE create_ewald_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: print_key, subsection

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="ewald", &
                          description="Ewald parameters controlling electrostatic only for CLASSICAL MM.", &
                          n_keywords=7, n_subsections=0, repeats=.FALSE., &
                          citations=(/Ewald1921, Darden1993, Essmann1995, Toukmaji1996, Laino2008/))

      NULLIFY (keyword, print_key, subsection)
      CALL keyword_create( &
         keyword, __LOCATION__, name="EWALD_TYPE", &
         description="The type of ewald you want to perform.", &
         citations=(/Ewald1921, Essmann1995, Darden1993/), &
         usage="EWALD_TYPE (NONE|EWALD|PME|SPME)", &
         default_i_val=do_ewald_ewald, &
         enum_c_vals=(/"none     ", &
                       "ewald    ", &
                       "pme      ", &
                       "spme     "/), &
         enum_i_vals=(/do_ewald_none, &
                       do_ewald_ewald, &
                       do_ewald_pme, &
                       do_ewald_spme/), &
         enum_desc=s2a("NONE standard real-space coulomb potential is computed together with the non-bonded contributions", &
                       "EWALD is the standard non-fft based ewald", &
                       "PME is the particle mesh using fft interpolation", &
                       "SPME is the smooth particle mesh using beta-Euler splines (recommended)"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="EWALD_ACCURACY", &
                          description="Expected accuracy in the Ewald sum. This number affects only the calculation of "// &
                          "the cutoff for the real-space term of the ewald summation (EWALD|PME|SPME) as well as the "// &
                          "construction of the neighbor lists (if the cutoff for non-bonded terms is smaller than the "// &
                          "value employed to compute the EWALD real-space term). This keyword has no "// &
                          "effect on the reciprocal space term (which can be tuned independently).", &
                          usage="EWALD_ACCURACY {real}", n_var=1, type_of_var=real_t, &
                          unit_str="hartree", default_r_val=1.0E-6_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="RCUT", &
                          description="Explicitly provide the real-space cutoff of the ewald summation (EWALD|PME|SPME). "// &
                          "This value is ignored in Tight-binding applications (rcut from basis overlap is used). "// &
                          "If present, overwrites the estimate of EWALD_ACCURACY and may affect the "// &
                          "construction of the neighbor lists for non-bonded terms (in FIST), if the value "// &
                          "specified is larger than the cutoff for non-bonded interactions.", &
                          usage="RCUT 5.0", n_var=1, type_of_var=real_t, unit_str="angstrom")
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="alpha", &
                          description="alpha parameter associated with Ewald (EWALD|PME|SPME). "// &
                          "Recommended for small systems is alpha = 3.5 / r_cut. "// &
                          "For Tight-binding application a recommended value is alpha = 1.0. "// &
                          "Tuning alpha, r_cut and gmax is needed to obtain O(N**1.5) scaling for ewald.", &
                          usage="alpha .30", &
                          default_r_val=cp_unit_to_cp2k(value=0.35_dp, unit_str="angstrom^-1"), &
                          unit_str='angstrom^-1')
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="gmax", &
                          description="number of grid points (SPME and EWALD). If a single number is specified, "// &
                          "the same number of points is used for all three directions on the grid. "// &
                          "If three numbers are given, each direction can have a different number of points. "// &
                          "The number of points needs to be FFTable (which depends on the library used) and odd for EWALD. "// &
                          "The optimal number depends e.g. on alpha and the size of the cell. 1 point per Angstrom is common.", &
                          usage="gmax 25 25 25", n_var=-1, type_of_var=integer_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="ns_max", &
                          description="number of grid points on small mesh (PME only), should be odd.", &
                          usage="ns_max 11", default_i_val=11)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="o_spline", &
                          description="order of the beta-Euler spline (SPME only)", &
                          usage="o_spline 6", default_i_val=6)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="epsilon", &
                          description="tolerance of gaussians for fft interpolation (PME only)", &
                          usage="epsilon 1e-6", default_r_val=1.e-6_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      NULLIFY (subsection)
      CALL create_rsgrid_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      NULLIFY (subsection)
      CALL section_create(subsection, __LOCATION__, name="MULTIPOLES", &
                          description="Enables the use of multipoles in the treatment of the electrostatics.", &
                          n_keywords=0, n_subsections=1, repeats=.FALSE., &
                          citations=(/Aguado2003, Laino2008/))

      CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
                          description="Controls the activation of the Multipoles", &
                          usage="&MULTIPOLES T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
      CALL section_add_keyword(subsection, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="MAX_MULTIPOLE_EXPANSION", &
                          description="Specify the maximum level of multipoles expansion used "// &
                          "for the electrostatics.", &
                          usage="MAX_MULTIPOLE_EXPANSION DIPOLE", &
                          enum_c_vals=s2a("NONE", "CHARGE", "DIPOLE", "QUADRUPOLE"), &
                          enum_desc=s2a("No multipolar terms! Check the codes providing a zero contribution.", &
                                        "Use up to the Charge term", &
                                        "Use up to the Dipole term", &
                                        "Use up to the Quadrupole term"), &
                          enum_i_vals=(/do_multipole_none, do_multipole_charge, do_multipole_dipole, &
                                        do_multipole_quadrupole/), type_of_var=enum_t)
      CALL section_add_keyword(subsection, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="POL_SCF", &
                          description="Specify the method to obtain self consistent induced "// &
                          "multipole moments.", &
                          usage="POL_SCF CONJUGATE_GRADIENT", &
                          enum_c_vals=s2a("NONE", "SELF_CONSISTENT", "CONJUGATE_GRADIENT"), &
                          enum_desc=s2a("No inducible multipoles.", &
                                        "Conventional self-consistent iteration.", &
                                        "Linear conjugate-gradient optimization of the sum "// &
                                        "of the electrostatic and induction energy. This "// &
                                        "method does not support non-linear polarization "// &
                                        "but is sometimes faster."), &
                          enum_i_vals=(/do_fist_pol_none, do_fist_pol_sc, do_fist_pol_cg/), &
                          type_of_var=enum_t, default_i_val=do_fist_pol_none)
      CALL section_add_keyword(subsection, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="MAX_IPOL_ITER", &
                          description="Specify the maximum number of iterations for induced "// &
                          "dipoles", &
                          usage="MAX_IPOL_ITER {int}", type_of_var=integer_t, &
                          n_var=1, default_i_val=0)
      CALL section_add_keyword(subsection, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="EPS_POL", &
                          description="Specify the rmsd threshold for the derivatives "// &
                          "of the energy towards the Cartesian dipoles components", &
                          usage="EPS_POL {real}", type_of_var=real_t, &
                          n_var=1, default_r_val=0.5e-07_dp)
      CALL section_add_keyword(subsection, keyword)
      CALL keyword_release(keyword)

      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      NULLIFY (subsection)
      CALL section_create(subsection, __LOCATION__, name="PRINT", &
                          description="Controls printing of Ewald properties", &
                          n_keywords=0, n_subsections=1, repeats=.FALSE.)
      NULLIFY (print_key)
      CALL cp_print_key_section_create(print_key, __LOCATION__, "PROGRAM_RUN_INFO", &
                                       description="controls the printing of ewald setup", &
                                       print_level=low_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
      CALL section_add_subsection(subsection, print_key)
      CALL section_release(print_key)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

   END SUBROUTINE create_ewald_section

! **************************************************************************************************
!> \brief creates the interpolation section for the periodic QM/MM
!> \param section ...
!> \author tlaino
! **************************************************************************************************
   SUBROUTINE create_gspace_interp_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: print_key

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="interpolator", &
                          description="controls the interpolation for the G-space term", &
                          n_keywords=5, n_subsections=0, repeats=.FALSE.)

      NULLIFY (keyword, print_key)

      CALL keyword_create(keyword, __LOCATION__, name="aint_precond", &
                          description="the approximate inverse to use to get the starting point"// &
                          " for the linear solver of the spline3 methods", &
                          usage="kind spline3", &
                          default_i_val=precond_spl3_aint, &
                          enum_c_vals=s2a("copy", "spl3_nopbc_aint1", "spl3_nopbc_precond1", &
                                          "spl3_nopbc_aint2", "spl3_nopbc_precond2", "spl3_nopbc_precond3"), &
                          enum_i_vals=(/no_precond, precond_spl3_aint, precond_spl3_1, &
                                        precond_spl3_aint2, precond_spl3_2, precond_spl3_3/))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="precond", &
                          description="The preconditioner used"// &
                          " for the linear solver of the spline3 methods", &
                          usage="kind spline3", &
                          default_i_val=precond_spl3_3, &
                          enum_c_vals=s2a("copy", "spl3_nopbc_aint1", "spl3_nopbc_precond1", &
                                          "spl3_nopbc_aint2", "spl3_nopbc_precond2", "spl3_nopbc_precond3"), &
                          enum_i_vals=(/no_precond, precond_spl3_aint, precond_spl3_1, &
                                        precond_spl3_aint2, precond_spl3_2, precond_spl3_3/))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="eps_x", &
                          description="accuracy on the solution for spline3 the interpolators", &
                          usage="eps_x 1.e-15", default_r_val=1.e-10_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="eps_r", &
                          description="accuracy on the residual for spline3 the interpolators", &
                          usage="eps_r 1.e-15", default_r_val=1.e-10_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="max_iter", &
                          variants=(/'maxiter'/), &
                          description="the maximum number of iterations", &
                          usage="max_iter 200", default_i_val=100)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      NULLIFY (print_key)
      CALL cp_print_key_section_create(print_key, __LOCATION__, "conv_info", &
                                       description="if convergence information about the linear solver"// &
                                       " of the spline methods should be printed", &
                                       print_level=medium_print_level, each_iter_names=s2a("SPLINE_FIND_COEFFS"), &
                                       each_iter_values=(/10/), filename="__STD_OUT__", &
                                       add_last=add_last_numeric)
      CALL section_add_subsection(section, print_key)
      CALL section_release(print_key)

   END SUBROUTINE create_gspace_interp_section

! **************************************************************************************************
!> \brief Creates the wavelet section
!> \param section the section to create
!> \author fschiff
!> \note
!>      this approach is based on the development of T. Deutsch and S. Goedecker
! **************************************************************************************************
   SUBROUTINE create_wavelet_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create( &
         section, __LOCATION__, name="wavelet", &
         description="Sets up parameters of wavelet based poisson solver.", &
         n_keywords=1, n_subsections=0, repeats=.FALSE., &
         citations=(/Genovese2006, Genovese2007/))

      NULLIFY (keyword)

      CALL keyword_create( &
         keyword, __LOCATION__, name="SCF_TYPE", &
         description="Type of scaling function used in the wavelet approach, the total energy depends on this choice, "// &
         "and the convergence with respect to cutoff depends on the selected scaling functions. "// &
         "Possible values are 8,14,16,20,24,30,40,50,60,100", &
         usage="SCF_TYPE integer", &
         n_var=1, default_i_val=40)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_wavelet_section

! **************************************************************************************************
!> \brief Creates the section for the implicit (generalized) poisson solver
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_implicit_ps_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: subsection

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="IMPLICIT", &
                          description="Parameters for the implicit (generalized) Poisson solver.", &
                          citations=(/BaniHashemian2016/), &
                          n_keywords=6, n_subsections=2, repeats=.FALSE.)

      NULLIFY (subsection, keyword)

      CALL create_dielectric_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_dbc_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL keyword_create( &
         keyword, __LOCATION__, name="BOUNDARY_CONDITIONS", &
         enum_c_vals=s2a('PERIODIC', 'MIXED', 'MIXED_PERIODIC', 'NEUMANN'), &
         enum_desc=s2a('periodic boundary conditions', 'Dirichlet + homogeneous Neumann boundary conditions', &
                       'Dirichlet + periodic boundary conditions', 'homogeneous Neumann BC (zero-average solution)'), &
         enum_i_vals=(/PERIODIC_BC, MIXED_BC, MIXED_PERIODIC_BC, NEUMANN_BC/), &
         description="Specifies the type of boundary conditions. Dirichlet=fixed value, Neumann=zero normal deriv. "// &
         "Mixed and Neumann boundaries essentially requires FFTW3 so that all grid sizes are FFT-able.", &
         usage="BOUNDARY_CONDITIONS <bc_type>", default_i_val=PERIODIC_BC)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="ZERO_INITIAL_GUESS", &
                          description="Whether or not to use zero potential as initial guess.", &
                          usage="ZERO_INITIAL_GUESS <logical>", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="max_iter", &
                          description="Maximum number of iterations.", &
                          usage="max_iter <integer>", default_i_val=30)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="tol", &
                          description="Stopping tolerance.", &
                          usage="tol <real>", default_r_val=1.0E-8_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="OR_PARAMETER", variants=s2a('omega'), &
                          description="Over-relaxation parameter (large epsilon requires smaller omega ~0.1).", &
                          usage="OR_PARAMETER <real>", default_r_val=1.0_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create( &
         keyword, __LOCATION__, name="NEUMANN_DIRECTIONS", &
         enum_c_vals=s2a('XYZ', 'XY', 'XZ', 'YZ', 'X', 'Y', 'Z'), &
         enum_i_vals=(/neumannXYZ, neumannXY, neumannXZ, neumannYZ, neumannX, neumannY, neumannZ/), &
         description="Directions in which homogeneous Neumann conditions are imposed. In the remaining directions "// &
         "periodic conditions will be enforced. Having specified MIXED or NEUMANN as BOUNDARY_CONDITIONS, "// &
         "the keyword is meant to be used to combine periodic and homogeneous Neumann conditions at the "// &
         "boundaries of the simulation cell.", &
         usage="NEUMANN_DIRECTIONS <direction>", default_i_val=neumannXYZ)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_implicit_ps_section

! **************************************************************************************************
!> \brief Creates the dielectric constant section.
!>  The dielectric constant is defined as a function of electronic density.
!>  [see O. Andreussi, I. Dabo, and N. Marzari, J. Chem. Phys., 136, 064102(2012)]
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_dielectric_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: subsection

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="DIELECTRIC", &
                          description="Parameters for the dielectric constant function.", &
                          n_keywords=6, n_subsections=2, repeats=.FALSE.)

      NULLIFY (keyword, subsection)

      CALL keyword_create(keyword, __LOCATION__, name="DIELECTRIC_CORE_CORRECTION", &
                          description="Avoid spurious values of the dielectric constant at the ionic core for pseudopotentials "// &
                          "where the electron density goes to zero at the core (e.g. GTH). "// &
                          "The correction is based on rho_core.", &
                          usage="DIELECTRIC_CORE_CORRECTION <logical>", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create( &
         keyword, __LOCATION__, name="DIELECTRIC_FUNCTION_TYPE", &
         enum_c_vals=s2a('density_dependent', 'spatially_dependent', 'spatially_rho_dependent'), &
         enum_i_vals=(/rho_dependent, spatially_dependent, spatially_rho_dependent/), &
         enum_desc=s2a("Dielectric constant as a function of the electron density "// &
                       "as e.g. proposed within the SCCS model.", &
                       "Various regions with different dielectric constants.", &
                       "Various regions with different dielectric constants. The dielectric constant decays to 1.0, "// &
                       "wherever the electron density is present."), &
         description="Preferred type for the dielectric constant function.", &
         usage="DIELECTRIC_FUNCTION_TYPE  <method>", default_i_val=rho_dependent)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="dielectric_constant", variants=s2a('epsilon'), &
                          description="Dielectric constant in the bulk of the solvent.", &
                          usage="dielectric_constant <real>", default_r_val=80.0_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="rho_min", &
                          description="Lower density threshold.", &
                          usage="rho_min <real>", default_r_val=1.0E-4_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="rho_max", &
                          description="Upper density threshold.", &
                          usage="rho_max <real>", default_r_val=1.0E-3_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create( &
         keyword, __LOCATION__, name="DERIVATIVE_METHOD", &
         enum_c_vals=s2a('fft', 'fft_use_deps', 'fft_use_drho', 'cd3', 'cd5', 'cd7'), &
         enum_i_vals=(/derivative_fft, derivative_fft_use_deps, derivative_fft_use_drho, &
                       derivative_cd3, derivative_cd5, derivative_cd7/), &
         enum_desc=s2a("FFT based deriv of epsilon, without correction (high cutoff needed).", &
                       "FFT based deriv of epsilon, with correction using gradient of epsilon (high cutoff needed).", &
                       "FFT based deriv of epsilon, with correction using gradient of rho (high cutoff needed).", &
                       "3-point central difference derivative.", &
                       "5-point central difference derivative.", &
                       "7-point central difference derivative (recommended)."), &
         description="Preferred method for evaluating the gradient of ln(eps).", &
         usage="DERIVATIVE_METHOD  <method>", default_i_val=derivative_cd7)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL create_dielec_aa_cuboidal_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_dielec_xaa_annular_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

   END SUBROUTINE create_dielectric_section

! **************************************************************************************************
!> \brief Creates the section for creating axis-aligned cuboidal dielectric region.
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_dielec_aa_cuboidal_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="DIELEC_AA_CUBOIDAL", &
                          description="Parameters for creating axis-aligned cuboidal dielectric region. "// &
                          "Note that once such a region is defined, the 'background' dielectric constant "// &
                          "would be the default (80.0), unless a different value is specified using the "// &
                          "keyword IMPLICIT%DIELECTRIC%DIELECTRIC_CONSTANT.", &
                          n_keywords=5, n_subsections=0, repeats=.TRUE.)

      NULLIFY (keyword)

      CALL keyword_create(keyword, __LOCATION__, name="dielectric_constant", variants=s2a('epsilon'), &
                          description="value of the dielectric constant inside the region.", &
                          usage="dielectric_constant <real>", default_r_val=80.0_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="X_xtnt", &
                          description="The X extents of the cuboid.", &
                          usage="X_xtnt <xmin(real)> <xmax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="Y_xtnt", &
                          description="The Y extents of the cuboid.", &
                          usage="Y_xtnt <ymin(real)> <ymax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="Z_xtnt", &
                          description="The Z extents of the cuboid.", &
                          usage="Z_xtnt <zmin(real)> <zmax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="SMOOTHING_WIDTH", variants=s2a('zeta'), &
                          description="The width of the standard mollifier.", &
                          usage="SMOOTHING_WIDTH <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.2_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_dielec_aa_cuboidal_section

! **************************************************************************************************
!> \brief Creates the section for creating x-axis-aligned annular dielectric region.
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_dielec_xaa_annular_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="DIELEC_XAA_ANNULAR", &
                          description="Parameters for creating x-axis-aligned annular dielectric region. "// &
                          "Note that once such a region is defined, the 'background' dielectric constant "// &
                          "would be the default (80.0), unless a different value is specified using the "// &
                          "keyword IMPLICIT%DIELECTRIC%DIELECTRIC_CONSTANT.", &
                          n_keywords=5, n_subsections=0, repeats=.TRUE.)

      NULLIFY (keyword)

      CALL keyword_create(keyword, __LOCATION__, name="dielectric_constant", variants=s2a('epsilon'), &
                          description="value of the dielectric constant inside the region.", &
                          usage="dielectric_constant <real>", default_r_val=80.0_dp)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="X_xtnt", &
                          description="The X extents of the annulus.", &
                          usage="X_xtnt <xmin(real)> <xmax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="base_center", &
                          description="The y and z coordinates of the annulus' base center.", &
                          usage="base_center <y(real)> <z(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="base_radii", &
                          description="The base radius of the annulus.", &
                          usage="base_radius <r1(real)> <r2(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="smoothing_width", variants=s2a('zeta'), &
                          description="The width of the standard mollifier.", &
                          usage="smoothing_width <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.2_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_dielec_xaa_annular_section

! **************************************************************************************************
!> \brief Creates the section for Dirichlet boundary conditions
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_dbc_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: subsection

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="DIRICHLET_BC", &
                          description="Parameters for creating Dirichlet type boundary conditions.", &
                          n_keywords=1, n_subsections=4, repeats=.FALSE.)

      NULLIFY (keyword)

      CALL keyword_create(keyword, __LOCATION__, name="VERBOSE_OUTPUT", &
                          description="Print out the coordinates of the vertices defining Dirichlet regions and their "// &
                          "tessellations (in Angstrom), the values of the electrostatic potential at the regions (in a.u.), "// &
                          "and their corresponding evaluated Lagrange multipliers.", &
                          usage="VERBOSE_OUTPUT <logical>", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      NULLIFY (subsection)

      CALL create_aa_planar_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_planar_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_aa_cylindrical_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

      CALL create_aa_cuboidal_section(subsection)
      CALL section_add_subsection(section, subsection)
      CALL section_release(subsection)

   END SUBROUTINE create_dbc_section

! **************************************************************************************************
!> \brief Creates the section for creating axis-aligned planar Dirichlet BC.
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_aa_planar_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="AA_PLANAR", &
                          description="Parameters for creating axis-aligned planar (rectangular) Dirichlet boundary regions.", &
                          n_keywords=10, n_subsections=0, repeats=.TRUE.)

      NULLIFY (keyword)

      CALL keyword_create(keyword, __LOCATION__, name="v_D", &
                          description="The value of the fixed potential to be imposed at the axis-aligned Dirichlet boundary.", &
                          usage="v_D <real>", unit_str="volt", type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="OSCILLATING_FRACTION", &
                          description="A fraction of the field can be set to oscilate over time.", &
                          usage="OSCILLATING_FRACTION <real>", default_r_val=0.0_dp, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="FREQUENCY", &
                          description="The frequency with which the oscillating fraction oscillates.", &
                          usage="FREQUENCY <real>", default_r_val=0.0_dp, unit_str="s^-1", type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PHASE", &
                          description="The phase of the oscillattion. A phase of zero corresponds to a cosine function. ", &
                          usage="PHASE <real>", default_r_val=0.0_dp, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PARALLEL_PLANE", &
                          enum_c_vals=s2a('XY', 'YZ', 'XZ'), &
                          enum_i_vals=(/xy_plane, yz_plane, xz_plane/), &
                          description="The coordinate plane that the region is parallel to.", &
                          usage="PARALLEL_PLANE <plane>", &
                          type_of_var=enum_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="INTERCEPT", &
                          description="The intercept of the rectangle's plane.", &
                          usage="INTERCEPT <real>", unit_str="angstrom", &
                          type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="X_xtnt", &
                          description="The X extents of the rectangle.", &
                          usage="X_xtnt <xmin(real)> <xmax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="Y_xtnt", &
                          description="The Y extents of the rectangle.", &
                          usage="Y_xtnt <ymin(real)> <ymax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="Z_xtnt", &
                          description="The Z extents of the rectangle.", &
                          usage="Z_xtnt <zmin(real)> <zmax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="N_PRTN", &
                          description="The number of partitions in the directions of the unit vectors generating the "// &
                          "corresponding PARALLEL_PLANE (e1, e2 or e3) for tiling the rectangluar region.", &
                          usage="N_PRTN <integer> <integer>", &
                          n_var=2, default_i_vals=(/1, 1/))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="THICKNESS", &
                          description="The thickness of the planar Dirichlet region.", &
                          usage="THICKNESS <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.75_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="SMOOTHING_WIDTH", variants=s2a('SIGMA'), &
                          description="The width of the transition region between the Dirichlet subdomain and its surrounding.", &
                          usage="SMOOTHING_WIDTH <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.2_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PERIODIC_REGION", &
                          description="Whether or not to take into consideration the effects of the periodicity of the "// &
                          "simulation cell (MIXED_PERIODIC bc) for regions defined sufficiently close to the boundaries.", &
                          usage="PERIODIC_REGION <logical>", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_aa_planar_section

! **************************************************************************************************
!> \brief Creates the section for creating axis-aligned planar Dirichlet BC.
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_planar_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create( &
         section, __LOCATION__, name="PLANAR", &
         description="Parameters for creating planar (rectangular) Dirichlet boundary regions with given vertices.", &
         n_keywords=7, n_subsections=0, repeats=.TRUE.)

      NULLIFY (keyword)

      CALL keyword_create(keyword, __LOCATION__, name="v_D", &
                          description="The value of the fixed potential to be imposed at the planar Dirichlet boundary.", &
                          usage="v_D <real>", unit_str="volt", type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="OSCILLATING_FRACTION", &
                          description="A fraction of the field can be set to oscilate over time.", &
                          usage="OSCILLATING_FRACTION <real>", default_r_val=0.0_dp, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="FREQUENCY", &
                          description="The frequency with which the oscillating fraction oscillates.", &
                          usage="FREQUENCY <real>", default_r_val=0.0_dp, unit_str="s^-1", type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PHASE", &
                          description="The phase of the oscillattion. A phase of zero corresponds to a cosine function. ", &
                          usage="PHASE <real>", default_r_val=0.0_dp, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="A", &
                          description="Coordinates of the vertex A.", &
                          usage="A <x(real)> <y(real)> <z(real)>", unit_str="angstrom", &
                          n_var=3, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="B", &
                          description="Coordinates of the vertex B.", &
                          usage="B <x(real)> <y(real)> <z(real)>", unit_str="angstrom", &
                          n_var=3, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="C", &
                          description="Coordinates of the vertex C.", &
                          usage="C <x(real)> <y(real)> <z(real)>", unit_str="angstrom", &
                          n_var=3, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create( &
         keyword, __LOCATION__, name="N_PRTN", &
         description="The number of partitions along the edges for tiling the rectangular region. If the edges "// &
         "have different lengths, from the two given values, the larger one will be assigned to the longer edge.", &
         usage="N_PRTN <integer> <integer>", &
         n_var=2, default_i_vals=(/1, 1/))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="THICKNESS", &
                          description="The thickness of the planar Dirichlet region.", &
                          usage="THICKNESS <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.75_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="SMOOTHING_WIDTH", variants=s2a('SIGMA'), &
                          description="The width of the transition region between the Dirichlet subdomain and its surrounding.", &
                          usage="SMOOTHING_WIDTH <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.2_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_planar_section

! **************************************************************************************************
!> \brief Creates the section for creating x-axis-aligned cylindrical Dirichlet BC.
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_aa_cylindrical_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="AA_CYLINDRICAL", &
                          description="Parameters for creating axis-aligned cylindrical Dirichlet boundary regions.", &
                          n_keywords=11, n_subsections=0, repeats=.TRUE.)

      NULLIFY (keyword)

      CALL keyword_create(keyword, __LOCATION__, name="v_D", &
                          description="The value of the fixed potential to be imposed at the cylindrical Dirichlet boundary.", &
                          usage="v_D <real>", unit_str="volt", type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="OSCILLATING_FRACTION", &
                          description="A fraction of the field can be set to oscilate over time.", &
                          usage="OSCILLATING_FRACTION <real>", default_r_val=0.0_dp, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="FREQUENCY", &
                          description="The frequency with which the oscillating fraction oscillates.", &
                          usage="FREQUENCY <real>", default_r_val=0.0_dp, unit_str="s^-1", type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PHASE", &
                          description="The phase of the oscillattion. A phase of zero corresponds to a cosine function. ", &
                          usage="PHASE <real>", default_r_val=0.0_dp, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PARALLEL_AXIS", &
                          enum_c_vals=s2a('X', 'Y', 'Z'), &
                          enum_i_vals=(/x_axis, y_axis, z_axis/), &
                          description="The coordinate axis that the cylindrical region extends along.", &
                          usage="PARALLEL_AXIS <axis>", &
                          type_of_var=enum_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="xtnt", &
                          description="The extents of the cylinder along its central axis.", &
                          usage="xtnt <xmin(real)> <xmax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="BASE_CENTER", &
                          description="The y and z coordinates (x and z or x and y coordinates, "// &
                          "depending on the choice of the parallel axis) of the cylinder's base center.", &
                          usage="BASE_CENTER <y(real)> <z(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="BASE_RADIUS", &
                          description="The base radius of the cylinder.", &
                          usage="BASE_RADIUS <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=1.0_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="N_SIDES", &
                          description="The number of sides (faces) of the n-gonal prism approximating the cylinder.", &
                          usage="N_SIDES <integer>", default_i_val=5)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="APX_TYPE", &
                          enum_c_vals=s2a('CIRCUMSCRIBED', 'INSCRIBED'), &
                          enum_i_vals=(/CIRCUMSCRIBED, INSCRIBED/), &
                          description="Specifies the type of the n-gonal prism approximating the cylinder.", &
                          usage="APX_TYPE <apx_type>", default_i_val=CIRCUMSCRIBED)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create( &
         keyword, __LOCATION__, name="N_PRTN", &
         description="The number of partitions along the face edges of the prism for tiling. If the edges "// &
         "have different lengths, from the two given values, the larger one will be assigned to the longer edge.", &
         usage="N_PRTN <integer> <integer>", &
         n_var=2, default_i_vals=(/1, 1/))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="THICKNESS", &
                          description="The thickness of the cylinder.", &
                          usage="THICKNESS <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.75_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="SMOOTHING_WIDTH", variants=s2a('SIGMA'), &
                          description="The width of the transition region between the Dirichlet subdomain and its surrounding.", &
                          usage="SMOOTHING_WIDTH <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.2_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create( &
         keyword, __LOCATION__, name="delta_alpha", &
         description="A central angle specifying the gap between the faces of the n-gonal prism. To avoide overlap "// &
         "between the cuboids (of the given thickness) built on top of the faces, a larger value is required if the"// &
         " number of faces (N_SIDES) is quite few and/or the base radius is fairly small.", &
         usage="delta_alpha <real>", default_r_val=0.05_dp, unit_str="rad")
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_aa_cylindrical_section

! **************************************************************************************************
!> \brief Creates the section for creating axis-aligned cuboidal Dirichlet region.
!> \param section the section to be created
!> \author Mohammad Hossein Bani-Hashemian
! **************************************************************************************************
   SUBROUTINE create_aa_cuboidal_section(section)
      TYPE(section_type), POINTER                        :: section

      TYPE(keyword_type), POINTER                        :: keyword

      CPASSERT(.NOT. ASSOCIATED(section))
      CALL section_create(section, __LOCATION__, name="AA_CUBOIDAL", &
                          description="Parameters for creating axis-aligned cuboidal Dirichlet regions.", &
                          n_keywords=7, n_subsections=0, repeats=.TRUE.)

      NULLIFY (keyword)

      CALL keyword_create(keyword, __LOCATION__, name="v_D", &
                          description="The value of the fixed potential to be imposed at the region.", &
                          usage="v_D <real>", unit_str="volt", type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="OSCILLATING_FRACTION", &
                          description="A fraction of the field can be set to oscilate over time.", &
                          usage="OSCILLATING_FRACTION <real>", default_r_val=0.0_dp, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="FREQUENCY", &
                          description="The frequency with which the oscillating fraction oscillates.", &
                          usage="FREQUENCY <real>", default_r_val=0.0_dp, unit_str="s^-1", type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PHASE", &
                          description="The phase of the oscillattion. A phase of zero corresponds to a cosine function. ", &
                          usage="PHASE <real>", default_r_val=0.0_dp, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="X_xtnt", &
                          description="The X extents of the cuboid.", &
                          usage="X_xtnt <xmin(real)> <xmax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="Y_xtnt", &
                          description="The Y extents of the cuboid.", &
                          usage="Y_xtnt <ymin(real)> <ymax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="Z_xtnt", &
                          description="The Z extents of the cuboid.", &
                          usage="Z_xtnt <zmin(real)> <zmax(real)>", unit_str="angstrom", &
                          n_var=2, type_of_var=real_t)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="N_PRTN", &
                          description="The number of partitions in the x, y and z directions for partitioning the cuboid.", &
                          usage="N_PRTN <integer> <integer> <integer>", &
                          n_var=3, default_i_vals=(/1, 1, 1/))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="SMOOTHING_WIDTH", variants=s2a('SIGMA'), &
                          description="The width of the transition region between the Dirichlet subdomain and its surrounding.", &
                          usage="SMOOTHING_WIDTH <real>", unit_str="angstrom", &
                          default_r_val=cp_unit_to_cp2k(value=0.2_dp, unit_str="angstrom"))
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

      CALL keyword_create(keyword, __LOCATION__, name="PERIODIC_REGION", &
                          description="Whether or not to take into consideration the effects of the periodicity of the "// &
                          "simulation cell (MIXED_PERIODIC bc) for regions defined sufficiently close to the boundaries.", &
                          usage="PERIODIC_REGION <logical>", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
      CALL section_add_keyword(section, keyword)
      CALL keyword_release(keyword)

   END SUBROUTINE create_aa_cuboidal_section

END MODULE input_cp2k_poisson
